Yampa
Brian
2/2/2023
NOTES FROM SF:
FOR STEAMBOAT, COMPARE THE DETRENDED CLIAMTE NORMAL CURVES AGAINST EACH OTHER
Derry & fassnacht 2010 Compare climatology of snotel statuions= they don’t have the same climatol9ogy based on the snow data
Figure 5b - https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2009WR007835
Disclaimer
This work is very preliminary as I get back into the coding swing of things. Data wrangling and figure generation will be done via R, but the rest of the project will be done using good ol’ microsoft products. This is just an entry point into data crunching and should by no means be considered a final product.
Also, I’m not great at this but whatever. I could automate this, but I’ll figure that out shortly!
Need to update figures & analyze sites for seasonal trends
Methodology
SNOTEL data was provided by the NRCS. Data was cleaned by removing outliers that are likely implausible; any year with more than 15 observations missing was removed. Temperatures were adjusted using the Morrisey method for stations identified by Ma et al (2019) due to SNOTEL temperature sensor changes, with the adjustment applied to pre-sensor change data. Daily mean observations were detrended to determine whether values were increasing or decreasing from the entire time series trend. Daily mean temperatures were first averaged by water year, with all water year means then averaged by day of water year. The mean temperature by day for the period of record was averaged. To find the standard deviation, the daily mean temperatures by water year was subtracted from the averaged mean temperature by day for the period of record. All water year means averaged by day of water year were subtracted from the temperature mean. The resulting values were then added together to find the “residual” of the daily mean temperatures by water year. The standard deviation was then computed from those residuals, with trends analyzed by Mann‐Kendall significance test and Theil‐Sen’s rate of change. Significant trends are identified with p-values of less than 0.10.
Morrisey Method
The Morrisey Method is taken from Ma, Fassnacht and Kampf..
In R script: T(adjusted) = 5.3x10(-7)xT(old)4+3.72x10(-5)xT(old)3-2.16x10(-3)xT(old)2-7.32x10^(-2)xT(old)+1.37
In the Ma et al. spreadsheet, H1 is Morrisey, H2 is Oiler
Yampa Area SNOTEL sites:
Burro Mountain 378 NSCE- Original, Bias- Morrisey 7/11/2002 (NSCE- 0.87 vs 0.84, Bias- -0.03 vs. 0.11)
Columbine 408 Morrisey 7/22/2005
Columbine Pass 409 Morrisey 6/23/2005
Crosho 426 Morrisey 7/21/2005
Dry Lake 457 Oyler -> Morrisey 7/30/2003
Elk River 467 Oyler -> Morrisey 8/7/2006
Lynx Pass 607 Morrisey 5/22/2006
Rabbit Ears 709 Morrisey 8/7/2006 (Does not go.)
Ripple Creek 717 Oyler -> Morrisey 8/7/2006
Tower 825 Original 8/18/2004
Trapper Lake 827 Original 12/13/2004
Data from thesis research:
SNOTEL_yampa_area <- snotel_download(site_id = c(378, 408, 409, 426, 457, 467, 607, 709, 717, 825, 827), path = tempdir('../data'), internal = TRUE)
write.csv(SNOTEL_yampa_area,"C:/Users/13074/Documents/ESS580/thesis_project/Yampa/data_raw/snotel_yampa.csv", row.names = FALSE) #write in the raw data
SNOTEL_yampa_area <- read.csv("C:/Users/13074/Documents/ESS580/thesis_project/Yampa/data_raw/snotel_yampa.csv", header = TRUE)Burro Mountain 378
NSCE- Original, Bias- Morrisey 7/11/2002 **Original is marked as having the smallest bias, but the sheet says original is smaller.)
snotel_378 <- SNOTEL_yampa_area %>%
filter(site_id == "378")#str(snotel_378) # check the date, usually a character.
snotel_378$Date <- as.Date(snotel_378$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_378_clean <- snotel_378 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_378_clean <- snotel_378_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))# Check for outliers
ggplot(snotel_378_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_378_clean <- snotel_378_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -40)# %>%
#filter(temperature_mean < 25)ggplot(snotel_378_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
#snotel_378_cull_count <- snotel_378_clean %>%
# filter(temperature_min > -40) %>%
# count(waterYear)
#snotel_378_cull_count
# filtering for too few observations in a year
snotel_378_cull_count_days <- snotel_378_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_378_cull_count_days## # A tibble: 10 x 2
## # Groups: waterYear [10]
## waterYear n
## <dbl> <int>
## 1 1985 1
## 2 1988 298
## 3 1989 341
## 4 1993 281
## 5 1994 261
## 6 1995 344
## 7 1996 320
## 8 2002 312
## 9 2016 310
## 10 2022 338snotel_378_clean_culled <- snotel_378_clean %>%
filter(waterYear != "1985" & waterYear != "1986" & waterYear != "1988" & waterYear != "1989" & waterYear != "1993" & waterYear != "1994" & waterYear != "1995" & waterYear != "1996" & waterYear != "2002" & waterYear != "2016" & waterYear != "2022")# & waterYear != "2017") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_378_clean_culled, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Culling WY 1986 as well.
ggplot(snotel_378_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_378_xts <- xts(snotel_378_clean_culled$temperature_mean, order.by = snotel_378_clean_culled$Date)
dygraph(temp_378_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") #snotel_378_clean_culled <- snotel_378_clean_culled %>%
# filter(temperature_mean < 19.3)
temp_378_xts <- xts(snotel_378_clean_culled$temperature_mean, order.by = snotel_378_clean_culled$Date)
dygraph(temp_378_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") Burro Mountain 378
NSCE- Original, Bias- Morrisey 7/11/2002
snotel_378_adjusted <- snotel_378_clean_culled %>%
mutate(temp_ad = if_else(Date < "2002-07-11", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))Non-corrected
378 Detrended
#using the clean culled df:
#average water year temperature
yearly_wy_aver_378 <- snotel_378_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_378 <- yearly_wy_aver_378 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_378 <- daily_wy_aver_378 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_378$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_378 <-daily_wy_aver_378 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_378$date_temp <- signif(daily_wy_aver2_378$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_378, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
378 SD
standard_dev_378 <- daily_wy_aver_378 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_378 <- standard_dev_378 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_378 <- standard_dev_all_378 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_378 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1987 | 8.562568 |
| 1990 | 9.114355 |
| 1991 | 9.254156 |
| 1992 | 8.275786 |
| 1997 | 8.898932 |
| 1998 | 8.983050 |
| 1999 | 8.073327 |
| 2000 | 8.686540 |
| 2001 | 9.301471 |
| 2003 | 8.509538 |
| 2004 | 8.017044 |
| 2005 | 7.723002 |
| 2006 | 8.547829 |
| 2007 | 8.796189 |
| 2008 | 8.749038 |
| 2009 | 7.771973 |
| 2010 | 8.617066 |
| 2011 | 8.282002 |
| 2012 | 8.445614 |
| 2013 | 8.991985 |
| 2014 | 8.166025 |
| 2015 | 7.657629 |
| 2017 | 8.034586 |
| 2018 | 8.020621 |
| 2019 | 8.518154 |
| 2020 | 8.848254 |
| 2021 | 8.718771 |
ggplot(standard_dev_all_378, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 378 average temperatures for water years 2005-2021
MK & SS for 378 (non-corrected)
sd_mk_378 <- mk.test(standard_dev_all_378$sd_2)
print(sd_mk_378)##
## Mann-Kendall trend test
##
## data: standard_dev_all_378$sd_2
## z = -1.5427, n = 27, p-value = 0.1229
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -75.0000000 2301.0000000 -0.2136752sd_sens_378 <- sens.slope(standard_dev_all_378$sd_2)
print(sd_sens_378)##
## Sen's slope
##
## data: standard_dev_all_378$sd_2
## z = -1.5427, n = 27, p-value = 0.1229
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.047456119 0.006621254
## sample estimates:
## Sen's slope
## -0.0230883Corrected
378 Detrended (corrected)
#using the clean culled df:
#average water year temperature
yearly_wy_aver_378_ad <- snotel_378_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_378_ad <- yearly_wy_aver_378_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_378_ad <- daily_wy_aver_378_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_378_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_378_ad <-daily_wy_aver_378_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_378_ad$date_temp_ad <- signif(daily_wy_aver2_378_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_378_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
378 SS (corrected)
standard_dev_378_ad <- daily_wy_aver_378_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_378_ad <- standard_dev_378_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_378_ad <- standard_dev_all_378_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_378_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1987 | 7.940547 |
| 1990 | 8.436056 |
| 1991 | 8.633602 |
| 1992 | 7.670874 |
| 1997 | 8.260208 |
| 1998 | 8.299486 |
| 1999 | 7.471734 |
| 2000 | 7.993287 |
| 2001 | 8.584377 |
| 2003 | 8.508522 |
| 2004 | 8.016141 |
| 2005 | 7.722598 |
| 2006 | 8.547499 |
| 2007 | 8.795362 |
| 2008 | 8.748133 |
| 2009 | 7.771479 |
| 2010 | 8.616338 |
| 2011 | 8.281602 |
| 2012 | 8.444783 |
| 2013 | 8.991116 |
| 2014 | 8.165236 |
| 2015 | 7.657576 |
| 2017 | 8.034186 |
| 2018 | 8.019986 |
| 2019 | 8.517552 |
| 2020 | 8.847387 |
| 2021 | 8.717958 |
ggplot(standard_dev_all_378_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 378 average temperatures for water years 1986-2021
MK & SS 378 (corrected)
sd_mk_378_ad <- mk.test(standard_dev_all_378_ad$sd_2)
print(sd_mk_378_ad)##
## Mann-Kendall trend test
##
## data: standard_dev_all_378_ad$sd_2
## z = 1.3342, n = 27, p-value = 0.1821
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 65.0000000 2301.0000000 0.1851852sd_sens_378_ad <- sens.slope(standard_dev_all_378_ad$sd_2)
print(sd_sens_378_ad)##
## Sen's slope
##
## data: standard_dev_all_378_ad$sd_2
## z = 1.3342, n = 27, p-value = 0.1821
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.009653315 0.035911433
## sample estimates:
## Sen's slope
## 0.01222034Columbine 408
Morrisey 7/22/2005
snotel_408 <- SNOTEL_yampa_area %>%
filter(site_id == "408")#str(snotel_408) # check the date, usually a character.
snotel_408$Date <- as.Date(snotel_408$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_408_clean <- snotel_408 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_408_clean <- snotel_408_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))# Check for outliers
ggplot(snotel_408_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_408_clean <- snotel_408_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -40)# %>%
#filter(temperature_mean < 25)ggplot(snotel_408_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
#snotel_408_cull_count <- snotel_408_clean %>%
# filter(temperature_min > -40) %>%
# count(waterYear)
#snotel_408_cull_count
# filtering for too few observations in a year
snotel_408_cull_count_days <- snotel_408_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_408_cull_count_days## # A tibble: 8 x 2
## # Groups: waterYear [8]
## waterYear n
## <dbl> <int>
## 1 1981 343
## 2 1983 317
## 3 1984 234
## 4 1985 195
## 5 1986 336
## 6 1987 298
## 7 1994 272
## 8 2002 347snotel_408_clean_culled <- snotel_408_clean %>%
filter(waterYear > "1987" & waterYear != "1993" & waterYear != "1994" & waterYear != "2002")# & waterYear != "1985" & waterYear != "1986" & waterYear != "1987" & waterYear != "1994" & waterYear != "2002")# & waterYear != "2002" & waterYear != "2016" & waterYear != "2022")# & waterYear != "2017") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_408_clean_culled, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Culling WY 1986 & 1993 as well.
ggplot(snotel_408_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_408_xts <- xts(snotel_408_clean_culled$temperature_mean, order.by = snotel_408_clean_culled$Date)
dygraph(temp_408_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") snotel_408_clean_culled <- snotel_408_clean_culled %>%
filter(temperature_mean > -30)
temp_408_xts <- xts(snotel_408_clean_culled$temperature_mean, order.by = snotel_408_clean_culled$Date)
dygraph(temp_408_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") Columbine 408 Morrisey 7/22/2005
snotel_408_adjusted <- snotel_408_clean_culled %>%
mutate(temp_ad = if_else(Date < "2005-07-22", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))Non-corrected
408 Detrended
#using the clean culled df:
#average water year temperature
yearly_wy_aver_408 <- snotel_408_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_408 <- yearly_wy_aver_408 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_408 <- daily_wy_aver_408 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_408$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_408 <-daily_wy_aver_408 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_408$date_temp <- signif(daily_wy_aver2_408$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_408, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
408 SD
standard_dev_408 <- daily_wy_aver_408 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_408 <- standard_dev_408 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_408 <- standard_dev_all_408 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_408 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1988 | 9.911311 |
| 1989 | 9.479817 |
| 1990 | 8.845588 |
| 1991 | 9.141713 |
| 1992 | 8.235717 |
| 1995 | 8.387811 |
| 1996 | 8.807923 |
| 1997 | 8.867338 |
| 1998 | 9.074387 |
| 1999 | 8.232888 |
| 2000 | 8.585091 |
| 2001 | 9.380546 |
| 2003 | 8.847298 |
| 2004 | 8.495495 |
| 2005 | 8.135796 |
| 2006 | 8.734304 |
| 2007 | 8.839901 |
| 2008 | 8.888992 |
| 2009 | 7.926154 |
| 2010 | 8.738675 |
| 2011 | 8.455577 |
| 2012 | 8.702054 |
| 2013 | 9.110505 |
| 2014 | 8.268271 |
| 2015 | 7.870630 |
| 2016 | 8.492743 |
| 2017 | 8.238399 |
| 2018 | 8.140089 |
| 2019 | 8.644103 |
| 2020 | 9.027143 |
| 2021 | 8.902353 |
| 2022 | 8.706784 |
ggplot(standard_dev_all_408, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 408 average temperatures for water years 2005-2021
MK & SS for 408 (non-corrected)
sd_mk_408 <- mk.test(standard_dev_all_408$sd_2)
print(sd_mk_408)##
## Mann-Kendall trend test
##
## data: standard_dev_all_408$sd_2
## z = -1.7676, n = 32, p-value = 0.07713
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -110.0000000 3802.6666667 -0.2217742sd_sens_408 <- sens.slope(standard_dev_all_408$sd_2)
print(sd_sens_408)##
## Sen's slope
##
## data: standard_dev_all_408$sd_2
## z = -1.7676, n = 32, p-value = 0.07713
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.033742130 0.001027723
## sample estimates:
## Sen's slope
## -0.01615257Corrected
408 Detrended (corrected)
#using the clean culled df:
#average water year temperature
yearly_wy_aver_408_ad <- snotel_408_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_408_ad <- yearly_wy_aver_408_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_408_ad <- daily_wy_aver_408_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_408_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_408_ad <-daily_wy_aver_408_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_408_ad$date_temp_ad <- signif(daily_wy_aver2_408_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_408_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
408 SS (corrected)
standard_dev_408_ad <- daily_wy_aver_408_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_408_ad <- standard_dev_408_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_408_ad <- standard_dev_all_408_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_408_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1988 | 9.255004 |
| 1989 | 8.873391 |
| 1990 | 8.216764 |
| 1991 | 8.545580 |
| 1992 | 7.664569 |
| 1995 | 7.780528 |
| 1996 | 8.186934 |
| 1997 | 8.260524 |
| 1998 | 8.410419 |
| 1999 | 7.651225 |
| 2000 | 7.928801 |
| 2001 | 8.691699 |
| 2003 | 8.174582 |
| 2004 | 7.895518 |
| 2005 | 7.482674 |
| 2006 | 8.735962 |
| 2007 | 8.841502 |
| 2008 | 8.890875 |
| 2009 | 7.928502 |
| 2010 | 8.741101 |
| 2011 | 8.457260 |
| 2012 | 8.703940 |
| 2013 | 9.112301 |
| 2014 | 8.270874 |
| 2015 | 7.872234 |
| 2016 | 8.495258 |
| 2017 | 8.239886 |
| 2018 | 8.142322 |
| 2019 | 8.646452 |
| 2020 | 9.028953 |
| 2021 | 8.904182 |
| 2022 | 8.708980 |
ggplot(standard_dev_all_408_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 408 average temperatures for water years 1986-2021
MK & SS 408 (corrected)
sd_mk_408_ad <- mk.test(standard_dev_all_408_ad$sd_2)
print(sd_mk_408_ad)##
## Mann-Kendall trend test
##
## data: standard_dev_all_408_ad$sd_2
## z = 1.2487, n = 32, p-value = 0.2118
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 78.0000000 3802.6666667 0.1572581sd_sens_408_ad <- sens.slope(standard_dev_all_408_ad$sd_2)
print(sd_sens_408_ad)##
## Sen's slope
##
## data: standard_dev_all_408_ad$sd_2
## z = 1.2487, n = 32, p-value = 0.2118
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.007349474 0.033195852
## sample estimates:
## Sen's slope
## 0.01316118Columbine Pass 409
Morrisey 6/23/2005
snotel_409 <- SNOTEL_yampa_area %>%
filter(site_id == "409")#str(snotel_409) # check the date, usually a character.
snotel_409$Date <- as.Date(snotel_409$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_409_clean <- snotel_409 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_409_clean <- snotel_409_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))# Check for outliers
ggplot(snotel_409_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_409_clean <- snotel_409_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -50) %>%
filter(temperature_mean < 40)ggplot(snotel_409_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
#snotel_409_cull_count <- snotel_409_clean %>%
# filter(temperature_min > -40) %>%
# count(waterYear)
#snotel_409_cull_count
# filtering for too few observations in a year
snotel_409_cull_count_days <- snotel_409_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_409_cull_count_days## # A tibble: 5 x 2
## # Groups: waterYear [5]
## waterYear n
## <dbl> <int>
## 1 1994 331
## 2 1995 190
## 3 2005 329
## 4 2019 348
## 5 2022 346snotel_409_clean_culled <- snotel_409_clean %>%
filter(waterYear != "1994" & waterYear != "1995" & waterYear != "2005" & waterYear != "2019" & waterYear != "2022")# & waterYear != "1986" & waterYear != "1987" & waterYear != "1994" & waterYear != "2002")# & waterYear != "2002" & waterYear != "2016" & waterYear != "2022")# & waterYear != "2017") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_409_clean_culled, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
ggplot(snotel_409_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_409_xts <- xts(snotel_409_clean_culled$temperature_mean, order.by = snotel_409_clean_culled$Date)
dygraph(temp_409_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") #snotel_409_clean_culled <- snotel_409_clean_culled %>%
# filter(temperature_mean > -30)
#temp_409_xts <- xts(snotel_409_clean_culled$temperature_mean, order.by = snotel_409_clean_culled$Date)
#dygraph(temp_409_xts) %>%
# dyAxis("y", label = "Daily mean temperature (°C)") Columbine Pass 409
Morrisey 6/23/2005
snotel_409_adjusted <- snotel_409_clean_culled %>%
mutate(temp_ad = if_else(Date < "2005-06-23", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))Non-corrected
409 Detrended
#using the clean culled df:
#average water year temperature
yearly_wy_aver_409 <- snotel_409_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_409 <- yearly_wy_aver_409 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_409 <- daily_wy_aver_409 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_409$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_409 <-daily_wy_aver_409 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_409$date_temp <- signif(daily_wy_aver2_409$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_409, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
409 SD
standard_dev_409 <- daily_wy_aver_409 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_409 <- standard_dev_409 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_409 <- standard_dev_all_409 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_409 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1987 | 8.326403 |
| 1988 | 9.167257 |
| 1989 | 8.983106 |
| 1990 | 8.555010 |
| 1991 | 7.839925 |
| 1992 | 7.987351 |
| 1993 | 8.471182 |
| 1996 | 8.508682 |
| 1997 | 8.478813 |
| 1998 | 8.929006 |
| 1999 | 7.623863 |
| 2000 | 8.631812 |
| 2001 | 9.047367 |
| 2002 | 9.516272 |
| 2003 | 8.791318 |
| 2004 | 8.708583 |
| 2006 | 8.125432 |
| 2007 | 8.536611 |
| 2008 | 8.816576 |
| 2009 | 7.846243 |
| 2010 | 8.511418 |
| 2011 | 8.496243 |
| 2012 | 8.285134 |
| 2013 | 8.792304 |
| 2014 | 7.958185 |
| 2015 | 7.393468 |
| 2016 | 8.187794 |
| 2017 | 8.066498 |
| 2018 | 7.805356 |
| 2020 | 8.856536 |
| 2021 | 8.674169 |
ggplot(standard_dev_all_409, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 409 average temperatures for water years 2005-2021
MK & SS for 409 (non-corrected)
sd_mk_409 <- mk.test(standard_dev_all_409$sd_2)
print(sd_mk_409)##
## Mann-Kendall trend test
##
## data: standard_dev_all_409$sd_2
## z = -1.1218, n = 31, p-value = 0.262
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -67.000000 3461.666667 -0.144086sd_sens_409 <- sens.slope(standard_dev_all_409$sd_2)
print(sd_sens_409)##
## Sen's slope
##
## data: standard_dev_all_409$sd_2
## z = -1.1218, n = 31, p-value = 0.262
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.034898600 0.009250748
## sample estimates:
## Sen's slope
## -0.01256069Corrected
409 Detrended (corrected)
#using the clean culled df:
#average water year temperature
yearly_wy_aver_409_ad <- snotel_409_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_409_ad <- yearly_wy_aver_409_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_409_ad <- daily_wy_aver_409_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_409_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_409_ad <-daily_wy_aver_409_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_409_ad$date_temp_ad <- signif(daily_wy_aver2_409_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_409_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
409 SS (corrected)
standard_dev_409_ad <- daily_wy_aver_409_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_409_ad <- standard_dev_409_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_409_ad <- standard_dev_all_409_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_409_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1987 | 7.660579 |
| 1988 | 8.456632 |
| 1989 | 8.306699 |
| 1990 | 7.862004 |
| 1991 | 7.209890 |
| 1992 | 7.340956 |
| 1993 | 7.806529 |
| 1996 | 7.809598 |
| 1997 | 7.810012 |
| 1998 | 8.200421 |
| 1999 | 7.004624 |
| 2000 | 7.910753 |
| 2001 | 8.302261 |
| 2002 | 8.749019 |
| 2003 | 8.046898 |
| 2004 | 8.036896 |
| 2006 | 8.125530 |
| 2007 | 8.536428 |
| 2008 | 8.816453 |
| 2009 | 7.846083 |
| 2010 | 8.511060 |
| 2011 | 8.495994 |
| 2012 | 8.285156 |
| 2013 | 8.792015 |
| 2014 | 7.958138 |
| 2015 | 7.393110 |
| 2016 | 8.187413 |
| 2017 | 8.066215 |
| 2018 | 7.805282 |
| 2020 | 8.856204 |
| 2021 | 8.673923 |
ggplot(standard_dev_all_409_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 409 average temperatures for water years 1986-2021
MK & SS 409 (corrected)
sd_mk_409_ad <- mk.test(standard_dev_all_409_ad$sd_2)
print(sd_mk_409_ad)##
## Mann-Kendall trend test
##
## data: standard_dev_all_409_ad$sd_2
## z = 2.1415, n = 31, p-value = 0.03223
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 127.0000000 3461.6666667 0.2731183sd_sens_409_ad <- sens.slope(standard_dev_all_409_ad$sd_2)
print(sd_sens_409_ad)##
## Sen's slope
##
## data: standard_dev_all_409_ad$sd_2
## z = 2.1415, n = 31, p-value = 0.03223
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## 0.002607708 0.039439754
## sample estimates:
## Sen's slope
## 0.02152587Crosho 426
Morrisey 7/21/2005
snotel_426 <- SNOTEL_yampa_area %>%
filter(site_id == "426")#str(snotel_426) # check the date, usually a character.
snotel_426$Date <- as.Date(snotel_426$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_426_clean <- snotel_426 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_426_clean <- snotel_426_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))# Check for outliers
ggplot(snotel_426_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_426_clean <- snotel_426_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -50) %>%
filter(temperature_mean < 40)ggplot(snotel_426_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
#snotel_426_cull_count <- snotel_426_clean %>%
# filter(temperature_min > -40) %>%
# count(waterYear)
#snotel_426_cull_count
# filtering for too few observations in a year
snotel_426_cull_count_days <- snotel_426_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_426_cull_count_days## # A tibble: 2 x 2
## # Groups: waterYear [2]
## waterYear n
## <dbl> <int>
## 1 2009 340
## 2 2021 349snotel_426_clean_culled <- snotel_426_clean %>%
filter(waterYear != "2009" & waterYear != "2021")# & waterYear != "2005" & waterYear != "2019" & waterYear != "2022")# & waterYear != "1986" & waterYear != "1987" & waterYear != "1994" & waterYear != "2002")# & waterYear != "2002" & waterYear != "2016" & waterYear != "2022")# & waterYear != "2017") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_426_clean_culled, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
ggplot(snotel_426_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_426_xts <- xts(snotel_426_clean_culled$temperature_mean, order.by = snotel_426_clean_culled$Date)
dygraph(temp_426_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") #snotel_426_clean_culled <- snotel_426_clean_culled %>%
# filter(temperature_mean > -30)
#temp_426_xts <- xts(snotel_426_clean_culled$temperature_mean, order.by = snotel_426_clean_culled$Date)
#dygraph(temp_426_xts) %>%
# dyAxis("y", label = "Daily mean temperature (°C)") Crosho 426 Morrisey 7/21/2005
snotel_426_adjusted <- snotel_426_clean_culled %>%
mutate(temp_ad = if_else(Date < "2005-07-21", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))Non-corrected
426 Detrended
#using the clean culled df:
#average water year temperature
yearly_wy_aver_426 <- snotel_426_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_426 <- yearly_wy_aver_426 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_426 <- daily_wy_aver_426 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_426$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_426 <-daily_wy_aver_426 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_426$date_temp <- signif(daily_wy_aver2_426$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_426, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
426 SD
standard_dev_426 <- daily_wy_aver_426 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_426 <- standard_dev_426 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_426 <- standard_dev_all_426 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_426 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1987 | 8.642717 |
| 1988 | 8.008775 |
| 1989 | 9.449036 |
| 1990 | 8.683053 |
| 1991 | 9.170577 |
| 1992 | 8.537791 |
| 1993 | 8.676514 |
| 1994 | 9.525407 |
| 1995 | 8.192424 |
| 1996 | 8.622041 |
| 1997 | 8.670823 |
| 1998 | 8.768147 |
| 1999 | 7.907919 |
| 2000 | 8.392275 |
| 2001 | 9.164446 |
| 2002 | 9.813834 |
| 2003 | 8.961799 |
| 2004 | 8.571169 |
| 2005 | 8.096051 |
| 2006 | 8.582421 |
| 2007 | 8.845773 |
| 2008 | 8.787455 |
| 2010 | 8.872551 |
| 2011 | 8.348895 |
| 2012 | 8.502421 |
| 2013 | 8.997153 |
| 2014 | 8.050838 |
| 2015 | 7.657484 |
| 2016 | 8.378217 |
| 2017 | 8.126881 |
| 2018 | 8.079436 |
| 2019 | 8.651915 |
| 2020 | 8.964956 |
| 2022 | 8.545430 |
ggplot(standard_dev_all_426, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 426 average temperatures for water years 2005-2021
MK & SS for 426 (non-corrected)
sd_mk_426 <- mk.test(standard_dev_all_426$sd_2)
print(sd_mk_426)##
## Mann-Kendall trend test
##
## data: standard_dev_all_426$sd_2
## z = -1.2749, n = 34, p-value = 0.2023
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -87.0000000 4550.3333333 -0.1550802sd_sens_426 <- sens.slope(standard_dev_all_426$sd_2)
print(sd_sens_426)##
## Sen's slope
##
## data: standard_dev_all_426$sd_2
## z = -1.2749, n = 34, p-value = 0.2023
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.02882424 0.00714949
## sample estimates:
## Sen's slope
## -0.009821212Corrected
426 Detrended (corrected)
#using the clean culled df:
#average water year temperature
yearly_wy_aver_426_ad <- snotel_426_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_426_ad <- yearly_wy_aver_426_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_426_ad <- daily_wy_aver_426_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_426_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_426_ad <-daily_wy_aver_426_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_426_ad$date_temp_ad <- signif(daily_wy_aver2_426_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_426_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
426 SS (corrected)
standard_dev_426_ad <- daily_wy_aver_426_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_426_ad <- standard_dev_426_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_426_ad <- standard_dev_all_426_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_426_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1987 | 8.032449 |
| 1988 | 7.469399 |
| 1989 | 8.830085 |
| 1990 | 8.043543 |
| 1991 | 8.571182 |
| 1992 | 7.909178 |
| 1993 | 8.068779 |
| 1994 | 8.814518 |
| 1995 | 7.580772 |
| 1996 | 7.973712 |
| 1997 | 8.044413 |
| 1998 | 8.100725 |
| 1999 | 7.332358 |
| 2000 | 7.718765 |
| 2001 | 8.465805 |
| 2002 | 9.079917 |
| 2003 | 8.257697 |
| 2004 | 7.940383 |
| 2005 | 7.444764 |
| 2006 | 8.582735 |
| 2007 | 8.846403 |
| 2008 | 8.788093 |
| 2010 | 8.872827 |
| 2011 | 8.349023 |
| 2012 | 8.502544 |
| 2013 | 8.997369 |
| 2014 | 8.051074 |
| 2015 | 7.658079 |
| 2016 | 8.378828 |
| 2017 | 8.126888 |
| 2018 | 8.080106 |
| 2019 | 8.652534 |
| 2020 | 8.965133 |
| 2022 | 8.546170 |
ggplot(standard_dev_all_426_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 426 average temperatures for water years 1986-2021
MK & SS 426 (corrected)
sd_mk_426_ad <- mk.test(standard_dev_all_426_ad$sd_2)
print(sd_mk_426_ad)##
## Mann-Kendall trend test
##
## data: standard_dev_all_426_ad$sd_2
## z = 1.9272, n = 34, p-value = 0.05396
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 131.0000000 4550.3333333 0.2335116sd_sens_426_ad <- sens.slope(standard_dev_all_426_ad$sd_2)
print(sd_sens_426_ad)##
## Sen's slope
##
## data: standard_dev_all_426_ad$sd_2
## z = 1.9272, n = 34, p-value = 0.05396
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.0008852415 0.0327236747
## sample estimates:
## Sen's slope
## 0.01613353Dry Lake 457
Oyler -> Morrisey 7/30/2003
snotel_457 <- SNOTEL_yampa_area %>%
filter(site_id == "457")#str(snotel_457) # check the date, usually a character.
snotel_457$Date <- as.Date(snotel_457$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_457_clean <- snotel_457 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_457_clean <- snotel_457_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))# Check for outliers
ggplot(snotel_457_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_457_clean <- snotel_457_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -40) %>%
filter(temperature_mean < 40)ggplot(snotel_457_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
#snotel_457_cull_count <- snotel_457_clean %>%
# filter(temperature_min > -40) %>%
# count(waterYear)
#snotel_457_cull_count
# filtering for too few observations in a year
snotel_457_cull_count_days <- snotel_457_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_457_cull_count_days## # A tibble: 6 x 2
## # Groups: waterYear [6]
## waterYear n
## <dbl> <int>
## 1 1985 1
## 2 1994 332
## 3 1996 349
## 4 1997 329
## 5 2003 346
## 6 2021 346snotel_457_clean_culled <- snotel_457_clean %>%
filter(waterYear != "1985" & waterYear != "1986" & waterYear != "1994" & waterYear != "1996" & waterYear != "1997" & waterYear != "2003" & waterYear != "2021")# & waterYear != "1987" & waterYear != "1994" & waterYear != "2002")# & waterYear != "2002" & waterYear != "2016" & waterYear != "2022")# & waterYear != "2017") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_457_clean_culled, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Also 1986
ggplot(snotel_457_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_457_xts <- xts(snotel_457_clean_culled$temperature_mean, order.by = snotel_457_clean_culled$Date)
dygraph(temp_457_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") #snotel_457_clean_culled <- snotel_457_clean_culled %>%
# filter(temperature_mean > -30)
#temp_457_xts <- xts(snotel_457_clean_culled$temperature_mean, order.by = snotel_457_clean_culled$Date)
#dygraph(temp_457_xts) %>%
# dyAxis("y", label = "Daily mean temperature (°C)") Dry Lake 457 Oyler -> Morrisey 7/30/2003
snotel_457_adjusted <- snotel_457_clean_culled %>%
mutate(temp_ad = if_else(Date < "2003-07-30", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))Non-corrected
457 Detrended
#using the clean culled df:
#average water year temperature
yearly_wy_aver_457 <- snotel_457_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_457 <- yearly_wy_aver_457 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_457 <- daily_wy_aver_457 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_457$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_457 <-daily_wy_aver_457 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_457$date_temp <- signif(daily_wy_aver2_457$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_457, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
457 SD
standard_dev_457 <- daily_wy_aver_457 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_457 <- standard_dev_457 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_457 <- standard_dev_all_457 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_457 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1987 | 8.892848 |
| 1988 | 9.953110 |
| 1989 | 9.757795 |
| 1990 | 9.176039 |
| 1991 | 9.421986 |
| 1992 | 8.507405 |
| 1993 | 8.878013 |
| 1995 | 8.358076 |
| 1998 | 9.150633 |
| 1999 | 8.440975 |
| 2000 | 8.736877 |
| 2001 | 9.645923 |
| 2002 | 10.040430 |
| 2004 | 8.051569 |
| 2005 | 7.863386 |
| 2006 | 8.768314 |
| 2007 | 9.036800 |
| 2008 | 8.900491 |
| 2009 | 7.976123 |
| 2010 | 8.677359 |
| 2011 | 8.610741 |
| 2012 | 8.805624 |
| 2013 | 9.212073 |
| 2014 | 8.350688 |
| 2015 | 7.879517 |
| 2016 | 8.589405 |
| 2017 | 8.326662 |
| 2018 | 8.323975 |
| 2019 | 8.654312 |
| 2020 | 9.071537 |
| 2022 | 8.555147 |
ggplot(standard_dev_all_457, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 457 average temperatures for water years 2005-2021
MK & SS for 457 (non-corrected)
sd_mk_457 <- mk.test(standard_dev_all_457$sd_2)
print(sd_mk_457)##
## Mann-Kendall trend test
##
## data: standard_dev_all_457$sd_2
## z = -2.3795, n = 31, p-value = 0.01734
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -141.0000000 3461.6666667 -0.3032258sd_sens_457 <- sens.slope(standard_dev_all_457$sd_2)
print(sd_sens_457)##
## Sen's slope
##
## data: standard_dev_all_457$sd_2
## z = -2.3795, n = 31, p-value = 0.01734
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.049787462 -0.005081117
## sample estimates:
## Sen's slope
## -0.02638274Corrected
457 Detrended (corrected)
#using the clean culled df:
#average water year temperature
yearly_wy_aver_457_ad <- snotel_457_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_457_ad <- yearly_wy_aver_457_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_457_ad <- daily_wy_aver_457_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_457_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_457_ad <-daily_wy_aver_457_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_457_ad$date_temp_ad <- signif(daily_wy_aver2_457_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_457_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
457 SS (corrected)
standard_dev_457_ad <- daily_wy_aver_457_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_457_ad <- standard_dev_457_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_457_ad <- standard_dev_all_457_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_457_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1987 | 8.200529 |
| 1988 | 9.184255 |
| 1989 | 9.050691 |
| 1990 | 8.442843 |
| 1991 | 8.739959 |
| 1992 | 7.843690 |
| 1993 | 8.207632 |
| 1995 | 7.712921 |
| 1998 | 8.438042 |
| 1999 | 7.811811 |
| 2000 | 8.022352 |
| 2001 | 8.894507 |
| 2002 | 9.283490 |
| 2004 | 8.051478 |
| 2005 | 7.863308 |
| 2006 | 8.768371 |
| 2007 | 9.036790 |
| 2008 | 8.900641 |
| 2009 | 7.976349 |
| 2010 | 8.677249 |
| 2011 | 8.611094 |
| 2012 | 8.805285 |
| 2013 | 9.211821 |
| 2014 | 8.350803 |
| 2015 | 7.879602 |
| 2016 | 8.589409 |
| 2017 | 8.326313 |
| 2018 | 8.323899 |
| 2019 | 8.654521 |
| 2020 | 9.071250 |
| 2022 | 8.555395 |
ggplot(standard_dev_all_457_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 457 average temperatures for water years 1986-2021
MK & SS 457 (corrected)
sd_mk_457_ad <- mk.test(standard_dev_all_457_ad$sd_2)
print(sd_mk_457_ad)##
## Mann-Kendall trend test
##
## data: standard_dev_all_457_ad$sd_2
## z = 0.44191, n = 31, p-value = 0.6586
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 2.700000e+01 3.461667e+03 5.806452e-02sd_sens_457_ad <- sens.slope(standard_dev_all_457_ad$sd_2)
print(sd_sens_457_ad)##
## Sen's slope
##
## data: standard_dev_all_457_ad$sd_2
## z = 0.44191, n = 31, p-value = 0.6586
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.01772587 0.02844930
## sample estimates:
## Sen's slope
## 0.004837867Elk River 467
Oyler -> Morrisey 8/7/2006
snotel_467 <- SNOTEL_yampa_area %>%
filter(site_id == "467")#str(snotel_467) # check the date, usually a character.
snotel_467$Date <- as.Date(snotel_467$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_467_clean <- snotel_467 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_467_clean <- snotel_467_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))# Check for outliers
ggplot(snotel_467_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_467_clean <- snotel_467_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -40) %>%
filter(temperature_mean < 40)ggplot(snotel_467_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
#snotel_467_cull_count <- snotel_467_clean %>%
# filter(temperature_min > -40) %>%
# count(waterYear)
#snotel_467_cull_count
# filtering for too few observations in a year
snotel_467_cull_count_days <- snotel_467_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_467_cull_count_days## # A tibble: 6 x 2
## # Groups: waterYear [6]
## waterYear n
## <dbl> <int>
## 1 1985 1
## 2 1994 326
## 3 1998 342
## 4 1999 281
## 5 2001 339
## 6 2021 344snotel_467_clean_culled <- snotel_467_clean %>%
filter(waterYear != "1985" & waterYear != "1986" & waterYear != "1994" & waterYear != "1998" & waterYear != "1999" & waterYear != "2001" & waterYear != "2021")# & waterYear != "2021")# & waterYear != "1987" & waterYear != "1994" & waterYear != "2002")# & waterYear != "2002" & waterYear != "2016" & waterYear != "2022")# & waterYear != "2017") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_467_clean_culled, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Also 1986
ggplot(snotel_467_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_467_xts <- xts(snotel_467_clean_culled$temperature_mean, order.by = snotel_467_clean_culled$Date)
dygraph(temp_467_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") #snotel_467_clean_culled <- snotel_467_clean_culled %>%
# filter(temperature_mean > -30)
#temp_467_xts <- xts(snotel_467_clean_culled$temperature_mean, order.by = snotel_467_clean_culled$Date)
#dygraph(temp_467_xts) %>%
# dyAxis("y", label = "Daily mean temperature (°C)") Elk River 467 Oyler -> Morrisey 8/7/2006
snotel_467_adjusted <- snotel_467_clean_culled %>%
mutate(temp_ad = if_else(Date < "2006-08-07", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))Non-corrected
467 Detrended
#using the clean culled df:
#average water year temperature
yearly_wy_aver_467 <- snotel_467_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_467 <- yearly_wy_aver_467 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_467 <- daily_wy_aver_467 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_467$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_467 <-daily_wy_aver_467 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_467$date_temp <- signif(daily_wy_aver2_467$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_467, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
467 SD
standard_dev_467 <- daily_wy_aver_467 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_467 <- standard_dev_467 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_467 <- standard_dev_all_467 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_467 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1987 | 9.103406 |
| 1988 | 10.510726 |
| 1989 | 10.049523 |
| 1990 | 9.547427 |
| 1991 | 9.692602 |
| 1992 | 8.717383 |
| 1993 | 9.114798 |
| 1995 | 8.625842 |
| 1996 | 8.985964 |
| 1997 | 9.119643 |
| 2000 | 8.850761 |
| 2002 | 10.120172 |
| 2003 | 9.161089 |
| 2004 | 8.589463 |
| 2005 | 8.525521 |
| 2006 | 9.564659 |
| 2007 | 9.088368 |
| 2008 | 9.090316 |
| 2009 | 8.102303 |
| 2010 | 8.886760 |
| 2011 | 8.734057 |
| 2012 | 8.902720 |
| 2013 | 9.375328 |
| 2014 | 8.477094 |
| 2015 | 8.044768 |
| 2016 | 8.660871 |
| 2017 | 8.514442 |
| 2018 | 8.465435 |
| 2019 | 8.883368 |
| 2020 | 9.345838 |
| 2022 | 8.902383 |
ggplot(standard_dev_all_467, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 467 average temperatures for water years 2005-2021
MK & SS for 467 (non-corrected)
sd_mk_467 <- mk.test(standard_dev_all_467$sd_2)
print(sd_mk_467)##
## Mann-Kendall trend test
##
## data: standard_dev_all_467$sd_2
## z = -2.7194, n = 31, p-value = 0.00654
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -161.0000000 3461.6666667 -0.3462366sd_sens_467 <- sens.slope(standard_dev_all_467$sd_2)
print(sd_sens_467)##
## Sen's slope
##
## data: standard_dev_all_467$sd_2
## z = -2.7194, n = 31, p-value = 0.00654
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.053188331 -0.007753413
## sample estimates:
## Sen's slope
## -0.02739628Corrected
467 Detrended (corrected)
#using the clean culled df:
#average water year temperature
yearly_wy_aver_467_ad <- snotel_467_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_467_ad <- yearly_wy_aver_467_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_467_ad <- daily_wy_aver_467_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_467_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_467_ad <-daily_wy_aver_467_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_467_ad$date_temp_ad <- signif(daily_wy_aver2_467_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_467_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
467 SS (corrected)
standard_dev_467_ad <- daily_wy_aver_467_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_467_ad <- standard_dev_467_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_467_ad <- standard_dev_all_467_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_467_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1987 | 8.416565 |
| 1988 | 9.721370 |
| 1989 | 9.348058 |
| 1990 | 8.816292 |
| 1991 | 9.016762 |
| 1992 | 8.057564 |
| 1993 | 8.438420 |
| 1995 | 7.963510 |
| 1996 | 8.304891 |
| 1997 | 8.454695 |
| 2000 | 8.126220 |
| 2002 | 9.357756 |
| 2003 | 8.423881 |
| 2004 | 7.952025 |
| 2005 | 7.847462 |
| 2006 | 8.843490 |
| 2007 | 9.088080 |
| 2008 | 9.090597 |
| 2009 | 8.102286 |
| 2010 | 8.886876 |
| 2011 | 8.734381 |
| 2012 | 8.902800 |
| 2013 | 9.375482 |
| 2014 | 8.477382 |
| 2015 | 8.045353 |
| 2016 | 8.661195 |
| 2017 | 8.514535 |
| 2018 | 8.465471 |
| 2019 | 8.883707 |
| 2020 | 9.345830 |
| 2022 | 8.902408 |
ggplot(standard_dev_all_467_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 467 average temperatures for water years 1986-2021
MK & SS 467 (corrected)
sd_mk_467_ad <- mk.test(standard_dev_all_467_ad$sd_2)
print(sd_mk_467_ad)##
## Mann-Kendall trend test
##
## data: standard_dev_all_467_ad$sd_2
## z = 0.57788, n = 31, p-value = 0.5633
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 3.500000e+01 3.461667e+03 7.526882e-02sd_sens_467_ad <- sens.slope(standard_dev_all_467_ad$sd_2)
print(sd_sens_467_ad)##
## Sen's slope
##
## data: standard_dev_all_467_ad$sd_2
## z = 0.57788, n = 31, p-value = 0.5633
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.01785969 0.02658485
## sample estimates:
## Sen's slope
## 0.003440178Lynx Pass 607
Morrisey 5/22/2006
snotel_607 <- SNOTEL_yampa_area %>%
filter(site_id == "607")#str(snotel_607) # check the date, usually a character.
snotel_607$Date <- as.Date(snotel_607$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_607_clean <- snotel_607 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_607_clean <- snotel_607_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))# Check for outliers
ggplot(snotel_607_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_607_clean <- snotel_607_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -40) %>%
filter(temperature_mean < 40)ggplot(snotel_607_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
#snotel_607_cull_count <- snotel_607_clean %>%
# filter(temperature_min > -40) %>%
# count(waterYear)
#snotel_607_cull_count
# filtering for too few observations in a year
snotel_607_cull_count_days <- snotel_607_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_607_cull_count_days## # A tibble: 4 x 2
## # Groups: waterYear [4]
## waterYear n
## <dbl> <int>
## 1 1985 1
## 2 1994 333
## 3 2019 284
## 4 2020 328snotel_607_clean_culled <- snotel_607_clean %>%
filter(waterYear != "1985" & waterYear != "1986" & waterYear != "1994" & waterYear != "2019" & waterYear != "2020")# & waterYear != "2001" & waterYear != "2021")# & waterYear != "2021")# & waterYear != "1987" & waterYear != "1994" & waterYear != "2002")# & waterYear != "2002" & waterYear != "2016" & waterYear != "2022")# & waterYear != "2017") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_607_clean_culled, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Also 1986
ggplot(snotel_607_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_607_xts <- xts(snotel_607_clean_culled$temperature_mean, order.by = snotel_607_clean_culled$Date)
dygraph(temp_607_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") #snotel_607_clean_culled <- snotel_607_clean_culled %>%
# filter(temperature_mean > -30)
#temp_607_xts <- xts(snotel_607_clean_culled$temperature_mean, order.by = snotel_607_clean_culled$Date)
#dygraph(temp_607_xts) %>%
# dyAxis("y", label = "Daily mean temperature (°C)") Lynx Pass 607 Morrisey 5/22/2006
snotel_607_adjusted <- snotel_607_clean_culled %>%
mutate(temp_ad = if_else(Date < "2006-05-22", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))Non-corrected
607 Detrended
#using the clean culled df:
#average water year temperature
yearly_wy_aver_607 <- snotel_607_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_607 <- yearly_wy_aver_607 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_607 <- daily_wy_aver_607 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_607$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_607 <-daily_wy_aver_607 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_607$date_temp <- signif(daily_wy_aver2_607$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_607, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
607 SD
standard_dev_607 <- daily_wy_aver_607 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_607 <- standard_dev_607 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_607 <- standard_dev_all_607 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_607 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1987 | 8.658921 |
| 1988 | 9.691007 |
| 1989 | 9.677093 |
| 1990 | 9.127534 |
| 1991 | 9.267600 |
| 1992 | 8.495617 |
| 1993 | 8.542784 |
| 1995 | 8.427463 |
| 1996 | 8.936524 |
| 1997 | 8.897439 |
| 1998 | 9.097517 |
| 1999 | 8.407095 |
| 2000 | 8.685669 |
| 2001 | 9.559462 |
| 2002 | 10.100422 |
| 2003 | 8.935295 |
| 2004 | 8.581025 |
| 2005 | 8.571462 |
| 2006 | 9.563477 |
| 2007 | 9.027055 |
| 2008 | 8.751128 |
| 2009 | 7.935985 |
| 2010 | 9.022507 |
| 2011 | 8.612122 |
| 2012 | 8.812590 |
| 2013 | 9.410619 |
| 2014 | 8.353556 |
| 2015 | 7.907493 |
| 2016 | 8.522596 |
| 2017 | 8.366193 |
| 2018 | 8.240463 |
| 2021 | 8.972977 |
| 2022 | 8.767136 |
ggplot(standard_dev_all_607, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 607 average temperatures for water years 2005-2021
MK & SS for 607 (non-corrected)
sd_mk_607 <- mk.test(standard_dev_all_607$sd_2)
print(sd_mk_607)##
## Mann-Kendall trend test
##
## data: standard_dev_all_607$sd_2
## z = -2.0298, n = 33, p-value = 0.04238
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -132.000 4165.333 -0.250sd_sens_607 <- sens.slope(standard_dev_all_607$sd_2)
print(sd_sens_607)##
## Sen's slope
##
## data: standard_dev_all_607$sd_2
## z = -2.0298, n = 33, p-value = 0.04238
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.036056261 -0.001515865
## sample estimates:
## Sen's slope
## -0.01788946Corrected
607 Detrended (corrected)
#using the clean culled df:
#average water year temperature
yearly_wy_aver_607_ad <- snotel_607_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_607_ad <- yearly_wy_aver_607_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_607_ad <- daily_wy_aver_607_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_607_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_607_ad <-daily_wy_aver_607_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_607_ad$date_temp_ad <- signif(daily_wy_aver2_607_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_607_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
607 SS (corrected)
standard_dev_607_ad <- daily_wy_aver_607_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_607_ad <- standard_dev_607_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_607_ad <- standard_dev_all_607_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_607_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1987 | 8.065805 |
| 1988 | 9.041423 |
| 1989 | 9.118850 |
| 1990 | 8.514101 |
| 1991 | 8.699806 |
| 1992 | 7.927056 |
| 1993 | 7.982533 |
| 1995 | 7.846909 |
| 1996 | 8.320426 |
| 1997 | 8.291259 |
| 1998 | 8.445612 |
| 1999 | 7.829298 |
| 2000 | 8.031514 |
| 2001 | 8.877490 |
| 2002 | 9.409575 |
| 2003 | 8.273822 |
| 2004 | 8.002022 |
| 2005 | 7.966244 |
| 2006 | 8.852307 |
| 2007 | 9.027252 |
| 2008 | 8.750916 |
| 2009 | 7.936108 |
| 2010 | 9.022503 |
| 2011 | 8.612197 |
| 2012 | 8.812993 |
| 2013 | 9.410822 |
| 2014 | 8.353805 |
| 2015 | 7.907862 |
| 2016 | 8.523025 |
| 2017 | 8.366000 |
| 2018 | 8.241006 |
| 2021 | 8.973005 |
| 2022 | 8.767834 |
ggplot(standard_dev_all_607_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 607 average temperatures for water years 1986-2021
MK & SS 607 (corrected)
sd_mk_607_ad <- mk.test(standard_dev_all_607_ad$sd_2)
print(sd_mk_607_ad)##
## Mann-Kendall trend test
##
## data: standard_dev_all_607_ad$sd_2
## z = 0.63527, n = 33, p-value = 0.5253
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 4.200000e+01 4.165333e+03 7.954545e-02sd_sens_607_ad <- sens.slope(standard_dev_all_607_ad$sd_2)
print(sd_sens_607_ad)##
## Sen's slope
##
## data: standard_dev_all_607_ad$sd_2
## z = 0.63527, n = 33, p-value = 0.5253
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.008825448 0.026678756
## sample estimates:
## Sen's slope
## 0.006477733Rabbit Ears 709
Morrisey 8/7/2006
snotel_709 <- SNOTEL_yampa_area %>%
filter(site_id == "709")#str(snotel_709) # check the date, usually a character.
snotel_709$Date <- as.Date(snotel_709$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_709_clean <- snotel_709 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_709_clean <- snotel_709_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))# Check for outliers
ggplot(snotel_709_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_709_clean <- snotel_709_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -40) %>%
filter(temperature_mean < 40)ggplot(snotel_709_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
#snotel_709_cull_count <- snotel_709_clean %>%
# filter(temperature_min > -40) %>%
# count(waterYear)
#snotel_709_cull_count
# filtering for too few observations in a year
snotel_709_cull_count_days <- snotel_709_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_709_cull_count_days## # A tibble: 14 x 2
## # Groups: waterYear [14]
## waterYear n
## <dbl> <int>
## 1 1992 318
## 2 1993 293
## 3 1994 332
## 4 1995 169
## 5 2006 348
## 6 2007 311
## 7 2011 292
## 8 2012 308
## 9 2013 300
## 10 2014 283
## 11 2015 44
## 12 2016 18
## 13 2017 34
## 14 2018 147snotel_709_clean_culled <- snotel_709_clean %>%
filter(waterYear > "1995" & waterYear != "2006" & waterYear != "2007" & waterYear != "2011" & waterYear != "2012" & waterYear != "2013" & waterYear != "2014" & waterYear != "2015" & waterYear != "2016" & waterYear != "2017" & waterYear != "2018")# & waterYear != "2002")# & waterYear != "2002" & waterYear != "2016" & waterYear != "2022")# & waterYear != "2017") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_709_clean_culled, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
ggplot(snotel_709_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_709_xts <- xts(snotel_709_clean_culled$temperature_mean, order.by = snotel_709_clean_culled$Date)
dygraph(temp_709_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") #snotel_709_clean_culled <- snotel_709_clean_culled %>%
# filter(temperature_mean > -30)
#temp_709_xts <- xts(snotel_709_clean_culled$temperature_mean, order.by = snotel_709_clean_culled$Date)
#dygraph(temp_709_xts) %>%
# dyAxis("y", label = "Daily mean temperature (°C)") Rabbit Ears 709 Morrisey 8/7/2006
snotel_709_adjusted <- snotel_709_clean_culled %>%
mutate(temp_ad = if_else(Date < "2006-08-07", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))Non-corrected
709 Detrended
#using the clean culled df:
#average water year temperature
yearly_wy_aver_709 <- snotel_709_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_709 <- yearly_wy_aver_709 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_709 <- daily_wy_aver_709 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_709$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_709 <-daily_wy_aver_709 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_709$date_temp <- signif(daily_wy_aver2_709$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_709, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
709 SD
standard_dev_709 <- daily_wy_aver_709 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_709 <- standard_dev_709 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_709 <- standard_dev_all_709 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_709 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1996 | 8.845436 |
| 1997 | 8.864608 |
| 1998 | 9.079136 |
| 1999 | 8.309577 |
| 2000 | 8.748503 |
| 2001 | 9.426260 |
| 2002 | 9.843776 |
| 2003 | 9.011443 |
| 2004 | 8.578795 |
| 2005 | 8.457080 |
| 2008 | 9.038776 |
| 2009 | 8.123905 |
| 2010 | 8.958458 |
| 2019 | 8.740431 |
| 2020 | 9.196584 |
| 2021 | 9.081772 |
| 2022 | 8.759943 |
ggplot(standard_dev_all_709, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 709 average temperatures for water years 2005-2021
MK & SS for 709 (non-corrected)
sd_mk_709 <- mk.test(standard_dev_all_709$sd_2)
print(sd_mk_709)##
## Mann-Kendall trend test
##
## data: standard_dev_all_709$sd_2
## z = 0.041193, n = 17, p-value = 0.9671
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 2.00000000 589.33333333 0.01470588sd_sens_709 <- sens.slope(standard_dev_all_709$sd_2)
print(sd_sens_709)##
## Sen's slope
##
## data: standard_dev_all_709$sd_2
## z = 0.041193, n = 17, p-value = 0.9671
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.04647212 0.03945218
## sample estimates:
## Sen's slope
## 0.0005780099Corrected
709 Detrended (corrected)
#using the clean culled df:
#average water year temperature
yearly_wy_aver_709_ad <- snotel_709_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_709_ad <- yearly_wy_aver_709_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_709_ad <- daily_wy_aver_709_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_709_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_709_ad <-daily_wy_aver_709_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_709_ad$date_temp_ad <- signif(daily_wy_aver2_709_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_709_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
709 SS (corrected)
standard_dev_709_ad <- daily_wy_aver_709_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_709_ad <- standard_dev_709_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_709_ad <- standard_dev_all_709_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_709_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1996 | 8.213076 |
| 1997 | 8.255603 |
| 1998 | 8.412792 |
| 1999 | 7.726856 |
| 2000 | 8.073670 |
| 2001 | 8.727586 |
| 2002 | 9.142221 |
| 2003 | 8.322566 |
| 2004 | 7.979937 |
| 2005 | 7.819602 |
| 2008 | 9.038980 |
| 2009 | 8.124255 |
| 2010 | 8.958543 |
| 2019 | 8.740727 |
| 2020 | 9.196244 |
| 2021 | 9.081146 |
| 2022 | 8.760212 |
ggplot(standard_dev_all_709_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 709 average temperatures for water years 1986-2021
MK & SS 709 (corrected)
sd_mk_709_ad <- mk.test(standard_dev_all_709_ad$sd_2)
print(sd_mk_709_ad)##
## Mann-Kendall trend test
##
## data: standard_dev_all_709_ad$sd_2
## z = 1.9361, n = 17, p-value = 0.05286
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 48.0000000 589.3333333 0.3529412sd_sens_709_ad <- sens.slope(standard_dev_all_709_ad$sd_2)
print(sd_sens_709_ad)##
## Sen's slope
##
## data: standard_dev_all_709_ad$sd_2
## z = 1.9361, n = 17, p-value = 0.05286
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.006786111 0.102902003
## sample estimates:
## Sen's slope
## 0.05101403Ripple Creek 717
Oyler -> Morrisey 8/7/2006
snotel_717 <- SNOTEL_yampa_area %>%
filter(site_id == "717")#str(snotel_717) # check the date, usually a character.
snotel_717$Date <- as.Date(snotel_717$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_717_clean <- snotel_717 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_717_clean <- snotel_717_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))# Check for outliers
ggplot(snotel_717_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_717_clean <- snotel_717_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -40) %>%
filter(temperature_mean < 40)ggplot(snotel_717_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
#snotel_717_cull_count <- snotel_717_clean %>%
# filter(temperature_min > -40) %>%
# count(waterYear)
#snotel_717_cull_count
# filtering for too few observations in a year
snotel_717_cull_count_days <- snotel_717_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_717_cull_count_days## # A tibble: 2 x 2
## # Groups: waterYear [2]
## waterYear n
## <dbl> <int>
## 1 1992 311
## 2 1993 330snotel_717_clean_culled <- snotel_717_clean %>%
filter(waterYear != "1992" & waterYear != "1993")# & waterYear != "2007" & waterYear != "2011" & waterYear != "2012" & waterYear != "2013" & waterYear != "2014" & waterYear != "2015" & waterYear != "2016" & waterYear != "2017" & waterYear != "2018")# & waterYear != "2002")# & waterYear != "2002" & waterYear != "2016" & waterYear != "2022")# & waterYear != "2017") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_717_clean_culled, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
ggplot(snotel_717_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_717_xts <- xts(snotel_717_clean_culled$temperature_mean, order.by = snotel_717_clean_culled$Date)
dygraph(temp_717_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") #snotel_717_clean_culled <- snotel_717_clean_culled %>%
# filter(temperature_mean > -30)
#temp_717_xts <- xts(snotel_717_clean_culled$temperature_mean, order.by = snotel_717_clean_culled$Date)
#dygraph(temp_717_xts) %>%
# dyAxis("y", label = "Daily mean temperature (°C)") Ripple Creek 717 Oyler -> Morrisey 8/7/2006
snotel_717_adjusted <- snotel_717_clean_culled %>%
mutate(temp_ad = if_else(Date < "2006-08-07", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))Non-corrected
717 Detrended
#using the clean culled df:
#average water year temperature
yearly_wy_aver_717 <- snotel_717_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_717 <- yearly_wy_aver_717 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_717 <- daily_wy_aver_717 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_717$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_717 <-daily_wy_aver_717 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_717$date_temp <- signif(daily_wy_aver2_717$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_717, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
717 SD
standard_dev_717 <- daily_wy_aver_717 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_717 <- standard_dev_717 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_717 <- standard_dev_all_717 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_717 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1987 | 8.803468 |
| 1988 | 9.707652 |
| 1989 | 9.270222 |
| 1990 | 9.163954 |
| 1991 | 9.043756 |
| 1994 | 9.475534 |
| 1995 | 8.216129 |
| 1996 | 8.766674 |
| 1997 | 8.761084 |
| 1998 | 8.941067 |
| 1999 | 8.106802 |
| 2000 | 8.669178 |
| 2001 | 9.209650 |
| 2002 | 9.676307 |
| 2003 | 9.008971 |
| 2004 | 8.362946 |
| 2005 | 8.373378 |
| 2006 | 9.248617 |
| 2007 | 8.800246 |
| 2008 | 8.921082 |
| 2009 | 7.998549 |
| 2010 | 8.681956 |
| 2011 | 8.605749 |
| 2012 | 8.582569 |
| 2013 | 9.046302 |
| 2014 | 8.329385 |
| 2015 | 7.737981 |
| 2016 | 8.402920 |
| 2017 | 8.250374 |
| 2018 | 8.202730 |
| 2019 | 8.544093 |
| 2020 | 8.992894 |
| 2021 | 8.957109 |
| 2022 | 8.519914 |
ggplot(standard_dev_all_717, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 717 average temperatures for water years 2005-2021
MK & SS for 717 (non-corrected)
sd_mk_717 <- mk.test(standard_dev_all_717$sd_2)
print(sd_mk_717)##
## Mann-Kendall trend test
##
## data: standard_dev_all_717$sd_2
## z = -2.6684, n = 34, p-value = 0.007621
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -181.0000000 4550.3333333 -0.3226381sd_sens_717 <- sens.slope(standard_dev_all_717$sd_2)
print(sd_sens_717)##
## Sen's slope
##
## data: standard_dev_all_717$sd_2
## z = -2.6684, n = 34, p-value = 0.007621
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.039224947 -0.006265457
## sample estimates:
## Sen's slope
## -0.02259117Corrected
717 Detrended (corrected)
#using the clean culled df:
#average water year temperature
yearly_wy_aver_717_ad <- snotel_717_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_717_ad <- yearly_wy_aver_717_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_717_ad <- daily_wy_aver_717_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_717_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_717_ad <-daily_wy_aver_717_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_717_ad$date_temp_ad <- signif(daily_wy_aver2_717_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_717_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
717 SS (corrected)
standard_dev_717_ad <- daily_wy_aver_717_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_717_ad <- standard_dev_717_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_717_ad <- standard_dev_all_717_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_717_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1987 | 8.203857 |
| 1988 | 9.063074 |
| 1989 | 8.683957 |
| 1990 | 8.529404 |
| 1991 | 8.477488 |
| 1994 | 8.832494 |
| 1995 | 7.643847 |
| 1996 | 8.162790 |
| 1997 | 8.181297 |
| 1998 | 8.315921 |
| 1999 | 7.566445 |
| 2000 | 8.038773 |
| 2001 | 8.563752 |
| 2002 | 9.013927 |
| 2003 | 8.355528 |
| 2004 | 7.809021 |
| 2005 | 7.779504 |
| 2006 | 8.600879 |
| 2007 | 8.799659 |
| 2008 | 8.920535 |
| 2009 | 7.997923 |
| 2010 | 8.681418 |
| 2011 | 8.605311 |
| 2012 | 8.581992 |
| 2013 | 9.045623 |
| 2014 | 8.328989 |
| 2015 | 7.737670 |
| 2016 | 8.402569 |
| 2017 | 8.249042 |
| 2018 | 8.202635 |
| 2019 | 8.543763 |
| 2020 | 8.992050 |
| 2021 | 8.956530 |
| 2022 | 8.520306 |
ggplot(standard_dev_all_717_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 717 average temperatures for water years 1986-2021
MK & SS 717 (corrected)
sd_mk_717_ad <- mk.test(standard_dev_all_717_ad$sd_2)
print(sd_mk_717_ad)##
## Mann-Kendall trend test
##
## data: standard_dev_all_717_ad$sd_2
## z = 0.56333, n = 34, p-value = 0.5732
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 3.900000e+01 4.550333e+03 6.951872e-02sd_sens_717_ad <- sens.slope(standard_dev_all_717_ad$sd_2)
print(sd_sens_717_ad)##
## Sen's slope
##
## data: standard_dev_all_717_ad$sd_2
## z = 0.56333, n = 34, p-value = 0.5732
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.01114957 0.02335425
## sample estimates:
## Sen's slope
## 0.004691148Tower 825
Original 8/18/2004
snotel_825 <- SNOTEL_yampa_area %>%
filter(site_id == "825")#str(snotel_825) # check the date, usually a character.
snotel_825$Date <- as.Date(snotel_825$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_825_clean <- snotel_825 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_825_clean <- snotel_825_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))# Check for outliers
ggplot(snotel_825_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_825_clean <- snotel_825_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -40) %>%
filter(temperature_mean < 30)ggplot(snotel_825_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
#snotel_825_cull_count <- snotel_825_clean %>%
# filter(temperature_min > -40) %>%
# count(waterYear)
#snotel_825_cull_count
# filtering for too few observations in a year
snotel_825_cull_count_days <- snotel_825_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_825_cull_count_days## # A tibble: 9 x 2
## # Groups: waterYear [9]
## waterYear n
## <dbl> <int>
## 1 1985 1
## 2 1987 348
## 3 1994 301
## 4 1995 347
## 5 1997 322
## 6 1998 343
## 7 1999 348
## 8 2002 327
## 9 2008 325snotel_825_clean_culled <- snotel_825_clean %>%
filter(waterYear > "1987" & waterYear != "1994" & waterYear != "1995" & waterYear != "1997" & waterYear != "1998" & waterYear != "1999" & waterYear != "2002" & waterYear != "2008")# & waterYear != "2017" & waterYear != "2018")# & waterYear != "2002")# & waterYear != "2002" & waterYear != "2016" & waterYear != "2022")# & waterYear != "2017") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_825_clean_culled, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
ggplot(snotel_825_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_825_xts <- xts(snotel_825_clean_culled$temperature_mean, order.by = snotel_825_clean_culled$Date)
dygraph(temp_825_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") #snotel_825_clean_culled <- snotel_825_clean_culled %>%
# filter(temperature_mean > -30)
#temp_825_xts <- xts(snotel_825_clean_culled$temperature_mean, order.by = snotel_825_clean_culled$Date)
#dygraph(temp_825_xts) %>%
# dyAxis("y", label = "Daily mean temperature (°C)") Tower 825 Original 8/18/2004
snotel_825_adjusted <- snotel_825_clean_culled %>%
mutate(temp_ad = if_else(Date < "2004-08-18", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))Non-corrected
825 Detrended
#using the clean culled df:
#average water year temperature
yearly_wy_aver_825 <- snotel_825_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_825 <- yearly_wy_aver_825 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_825 <- daily_wy_aver_825 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_825$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_825 <-daily_wy_aver_825 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_825$date_temp <- signif(daily_wy_aver2_825$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_825, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
825 SD
standard_dev_825 <- daily_wy_aver_825 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_825 <- standard_dev_825 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_825 <- standard_dev_all_825 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_825 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1988 | 10.158111 |
| 1989 | 9.521283 |
| 1990 | 9.387663 |
| 1991 | 9.495098 |
| 1992 | 8.468342 |
| 1993 | 8.686860 |
| 1996 | 9.236157 |
| 2000 | 8.951568 |
| 2001 | 9.626977 |
| 2003 | 9.317887 |
| 2004 | 8.794074 |
| 2005 | 8.168991 |
| 2006 | 9.046952 |
| 2007 | 9.112464 |
| 2009 | 8.322586 |
| 2010 | 8.961897 |
| 2011 | 8.871929 |
| 2012 | 8.997194 |
| 2013 | 9.382743 |
| 2014 | 8.649080 |
| 2015 | 8.104451 |
| 2016 | 8.792359 |
| 2017 | 8.593127 |
| 2018 | 8.634239 |
| 2019 | 8.862025 |
| 2020 | 9.328094 |
| 2021 | 9.270956 |
| 2022 | 9.124264 |
ggplot(standard_dev_all_825, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 825 average temperatures for water years 2005-2021
MK & SS for 825 (non-corrected)
sd_mk_825 <- mk.test(standard_dev_all_825$sd_2)
print(sd_mk_825)##
## Mann-Kendall trend test
##
## data: standard_dev_all_825$sd_2
## z = -1.7188, n = 28, p-value = 0.08565
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -88.0000000 2562.0000000 -0.2328042sd_sens_825 <- sens.slope(standard_dev_all_825$sd_2)
print(sd_sens_825)##
## Sen's slope
##
## data: standard_dev_all_825$sd_2
## z = -1.7188, n = 28, p-value = 0.08565
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.04047480 0.00485368
## sample estimates:
## Sen's slope
## -0.02125455Corrected
825 Detrended (corrected)
#using the clean culled df:
#average water year temperature
yearly_wy_aver_825_ad <- snotel_825_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_825_ad <- yearly_wy_aver_825_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_825_ad <- daily_wy_aver_825_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_825_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_825_ad <-daily_wy_aver_825_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_825_ad$date_temp_ad <- signif(daily_wy_aver2_825_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_825_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
825 SS (corrected)
standard_dev_825_ad <- daily_wy_aver_825_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_825_ad <- standard_dev_825_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_825_ad <- standard_dev_all_825_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_825_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1988 | 9.527196 |
| 1989 | 8.967107 |
| 1990 | 8.793798 |
| 1991 | 8.946793 |
| 1992 | 7.952054 |
| 1993 | 8.156645 |
| 1996 | 8.650702 |
| 2000 | 8.338506 |
| 2001 | 8.985410 |
| 2003 | 8.675745 |
| 2004 | 8.209448 |
| 2005 | 8.172186 |
| 2006 | 9.049660 |
| 2007 | 9.115791 |
| 2009 | 8.324666 |
| 2010 | 8.964371 |
| 2011 | 8.875058 |
| 2012 | 9.001462 |
| 2013 | 9.386059 |
| 2014 | 8.652436 |
| 2015 | 8.107386 |
| 2016 | 8.796286 |
| 2017 | 8.596142 |
| 2018 | 8.637340 |
| 2019 | 8.865259 |
| 2020 | 9.331765 |
| 2021 | 9.273785 |
| 2022 | 9.127176 |
ggplot(standard_dev_all_825_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 825 average temperatures for water years 1986-2021
MK & SS 825 (corrected)
sd_mk_825_ad <- mk.test(standard_dev_all_825_ad$sd_2)
print(sd_mk_825_ad)##
## Mann-Kendall trend test
##
## data: standard_dev_all_825_ad$sd_2
## z = 1.0076, n = 28, p-value = 0.3137
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 52.0000000 2562.0000000 0.1375661sd_sens_825_ad <- sens.slope(standard_dev_all_825_ad$sd_2)
print(sd_sens_825_ad)##
## Sen's slope
##
## data: standard_dev_all_825_ad$sd_2
## z = 1.0076, n = 28, p-value = 0.3137
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.009639429 0.034852064
## sample estimates:
## Sen's slope
## 0.01232876Trapper Lake 827 Original 12/13/2004
Trapper Lake 827
Original 12/13/2004
snotel_827 <- SNOTEL_yampa_area %>%
filter(site_id == "827")#str(snotel_827) # check the date, usually a character.
snotel_827$Date <- as.Date(snotel_827$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_827_clean <- snotel_827 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_827_clean <- snotel_827_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))# Check for outliers
ggplot(snotel_827_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_827_clean <- snotel_827_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -40) %>%
filter(temperature_mean < 30)ggplot(snotel_827_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
#snotel_827_cull_count <- snotel_827_clean %>%
# filter(temperature_min > -40) %>%
# count(waterYear)
#snotel_827_cull_count
# filtering for too few observations in a year
snotel_827_cull_count_days <- snotel_827_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_827_cull_count_days## # A tibble: 18 x 2
## # Groups: waterYear [18]
## waterYear n
## <dbl> <int>
## 1 1986 331
## 2 1994 179
## 3 1995 323
## 4 1996 53
## 5 1997 278
## 6 1999 341
## 7 2001 348
## 8 2002 310
## 9 2003 340
## 10 2004 342
## 11 2009 342
## 12 2010 346
## 13 2011 326
## 14 2013 236
## 15 2014 241
## 16 2016 321
## 17 2017 219
## 18 2019 332snotel_827_clean_culled <- snotel_827_clean %>%
filter(waterYear != "1986" & waterYear != "1994" & waterYear != "1995" & waterYear != "1996" & waterYear != "1997" & waterYear != "1999" & waterYear != "2001" & waterYear != "2002" & waterYear != "2003" & waterYear != "2004" & waterYear != "2009" & waterYear != "2010" & waterYear != "2011" & waterYear != "2013" & waterYear != "2014" & waterYear != "2016" & waterYear != "2017" & waterYear != "2019") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_827_clean_culled, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
ggplot(snotel_827_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_827_xts <- xts(snotel_827_clean_culled$temperature_mean, order.by = snotel_827_clean_culled$Date)
dygraph(temp_827_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") #snotel_827_clean_culled <- snotel_827_clean_culled %>%
# filter(temperature_mean > -30)
#temp_827_xts <- xts(snotel_827_clean_culled$temperature_mean, order.by = snotel_827_clean_culled$Date)
#dygraph(temp_827_xts) %>%
# dyAxis("y", label = "Daily mean temperature (°C)") Trapper Lake 827
Original 12/13/2004
snotel_827_adjusted <- snotel_827_clean_culled %>%
mutate(temp_ad = if_else(Date < "2004-12-13", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))Non-corrected
827 Detrended
#using the clean culled df:
#average water year temperature
yearly_wy_aver_827 <- snotel_827_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_827 <- yearly_wy_aver_827 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_827 <- daily_wy_aver_827 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_827$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_827 <-daily_wy_aver_827 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_827$date_temp <- signif(daily_wy_aver2_827$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_827, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
827 SD
standard_dev_827 <- daily_wy_aver_827 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_827 <- standard_dev_827 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_827 <- standard_dev_all_827 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_827 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1987 | 8.785352 |
| 1988 | 9.643636 |
| 1989 | 9.399426 |
| 1990 | 9.017517 |
| 1991 | 9.053817 |
| 1992 | 8.461780 |
| 1993 | 8.592203 |
| 1998 | 9.226580 |
| 2000 | 8.675610 |
| 2005 | 7.901837 |
| 2006 | 8.811208 |
| 2007 | 8.956871 |
| 2008 | 9.053611 |
| 2012 | 8.641372 |
| 2015 | 7.855779 |
| 2018 | 8.196461 |
| 2020 | 8.933495 |
| 2021 | 8.790312 |
| 2022 | 8.662496 |
ggplot(standard_dev_all_827, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 827 average temperatures for water years 2005-2021
MK & SS for 827 (non-corrected)
sd_mk_827 <- mk.test(standard_dev_all_827$sd_2)
print(sd_mk_827)##
## Mann-Kendall trend test
##
## data: standard_dev_all_827$sd_2
## z = -1.7493, n = 19, p-value = 0.08024
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -51.0000000 817.0000000 -0.2982456sd_sens_827 <- sens.slope(standard_dev_all_827$sd_2)
print(sd_sens_827)##
## Sen's slope
##
## data: standard_dev_all_827$sd_2
## z = -1.7493, n = 19, p-value = 0.08024
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.06891402 0.00422489
## sample estimates:
## Sen's slope
## -0.03456914Corrected
827 Detrended (corrected)
#using the clean culled df:
#average water year temperature
yearly_wy_aver_827_ad <- snotel_827_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_827_ad <- yearly_wy_aver_827_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_827_ad <- daily_wy_aver_827_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_827_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_827_ad <-daily_wy_aver_827_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_827_ad$date_temp_ad <- signif(daily_wy_aver2_827_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_827_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
827 SS (corrected)
standard_dev_827_ad <- daily_wy_aver_827_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_827_ad <- standard_dev_827_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_827_ad <- standard_dev_all_827_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_827_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1987 | 8.205336 |
| 1988 | 9.015708 |
| 1989 | 8.826659 |
| 1990 | 8.399945 |
| 1991 | 8.491804 |
| 1992 | 7.896678 |
| 1993 | 8.029814 |
| 1998 | 8.583994 |
| 2000 | 8.036502 |
| 2005 | 7.717045 |
| 2006 | 8.811366 |
| 2007 | 8.957047 |
| 2008 | 9.053812 |
| 2012 | 8.641373 |
| 2015 | 7.856092 |
| 2018 | 8.196931 |
| 2020 | 8.933512 |
| 2021 | 8.790758 |
| 2022 | 8.663052 |
ggplot(standard_dev_all_827_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 827 average temperatures for water years 1986-2021
MK & SS 827 (corrected)
sd_mk_827_ad <- mk.test(standard_dev_all_827_ad$sd_2)
print(sd_mk_827_ad)##
## Mann-Kendall trend test
##
## data: standard_dev_all_827_ad$sd_2
## z = 0.41983, n = 19, p-value = 0.6746
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 13.00000000 817.00000000 0.07602339sd_sens_827_ad <- sens.slope(standard_dev_all_827_ad$sd_2)
print(sd_sens_827_ad)##
## Sen's slope
##
## data: standard_dev_all_827_ad$sd_2
## z = 0.41983, n = 19, p-value = 0.6746
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.03007505 0.05895184
## sample estimates:
## Sen's slope
## 0.012232Note: figure captions are incorrect.